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Introduction 





The Large Hadron Collider is now at its third year of successful operation and both ATLAS 
and CMS report tantalizing hints of a Higgs boson at about 125 GeV. By the end of the 
2012 run the experiments are likely to be able to either confirm those hints as a firm 
discovery or else exclude any Standard Model (SM) Higgs boson. In the event of a firm 
discovery further detailed examination of various production and decay channels will be 
necessary to determine the nature of the Higgs sector. 

The dominant production channel in the SM, but also in all non-fermiophobic models of 
new physics, is single Higgs hadroproduction. Within the SM the production mechanism is 
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dominated by gluon fusion, since the alternative mechanism of quark annihilation is severely 
suppressed by the small Yukawa coupling of bottom and light quarks to the Higgs boson. 
However, if the Higgs sector is non-minimal, as is the case in any two-Higgs-doublet model 
(among which the MSSM is the most studied example) the Yukawa coupling to down-type 
quarks is enhanced by a factor of tan f3 (the ratio of the vacuum expectation values of the 
two doublets) and the contribution of the bb ^ H process becomes significant. Furthermore 
the production cross section through gluon fusion decreases due to the enhanced, negative 
top-bottom interference diagrams. In such a scenario, the production of a Higgs boson via 
bb pairs contributes much more than in the SM, and a detailed description of this process 
is desirable. In other BSM models, for example in models with dynamically generated 
Yukawa couplings [1, 2], both the bottom and the charm quarks have enhanced couplings 
to the Higgs boson and charm annihilation becomes important as well. 

The experimental searches are currently focused on measuring an enhanced production 
rate via bottom annihilation in the t~^t~ decay channel with the MSSM as the default 
BSM model [3, 4]. There are, moreover, several studies on measuring single Higgs decaying 
to bottom quarks in more generic models in which bottom annihilation is the dominant 
production channel, using, for example, three 6-tagged jets [5, 6], or measuring the ratio 
of three heavy (c- or b-) jet events to three b-jet events to discriminate between classes of 
models with two Higgs doublets [7]. 

Bottom quark annihilation has been the subject of much theoretical discussion in the 
last decade due to the freedom in treating the initial state bottom quarks. Bottom quarks lie 
in an intermediate mass range between the non-perturbative regime of the proton mass and 
the typical scale of a hard scattering event at the LHC. One can retain their small mass in 
the calculation, and exclude them from the proton constituents (four flavor scheme - 4FS) 
or treat them as massless partons with their own parton distribution functions (five flavor 
scheme - 5FS). In the 4FS the inclusive cross section develops large logarithms ~ log(^^) 
due to the collinear production of 6-quarks which is regulated by the bottom mass. In 
the 5FS these logarithms are re-summed to all orders by the DGLAP evolution inside 
the bottom PDFs, for all scales up to the factorization scale adopted in the calculation. 
Improved convergence of the perturbative expansion is an advantage of the 5FS approach, 
but at the same time it makes the 5FS prediction very sensitive to the choice of factorization 
scale. It has been realized that if the factorization scale is set to low values ~ m/f/4, both 
the 5FS and the 4FS predictions for the inclusive cross sections agree with each other 
within their respective uncertainties [8-10], and there is an open discussion as to how one 
would combine information from both approaches [11, 12]. 

In the 4FS, the lowest order process would be gg — t- bbH which begins at order in the 
QCD perturbative expansion and is known at next-to-leading-order (NLO) in QCD [13-16]. 
The process bg — >• bH, which starts at order q^, has also been studied at NLO in QCD [17] 
and with electroweak (EW) corrections [18]. In the 5FS the lowest process is bb — t- H. 
Hence the LO 4FS process where a non-collinear bottom pair is observable, is only reached 
for the first time at NNLO in the 5FS. The inclusive cross section, in the 5FS, of 66 — )• is 
known at NNLO in QCD [19] as well as at NLO in EW [20]. NNNLO threshold re-summed 
soft and collinear terms are also known [21] and the transverse momentum distribution 
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of the Higgs boson has been studied with re-summation techniques [22, 23] . Also known 
at NNLO are the zero-, one- and two-jet rates and related distributions [24], quantities 
which can be obtained already from the differential H + jet computation at NLO [25] in 
combination with the fully inclusive NNLO cross section. 

In this paper we present the fully differential NNLO QCD cross section for hb ^ H 
in the 5FS within the SM. NLO computations are currently performed with very well 
automated methods. Obtaining fully differential cross sections and decay rates at one order 
higher in the perturbative expansion requires the solution of new challenging problems. 
Regarding the treatment of the real emissions, pioneered for NLO computations in [26, 27], 
rapid progress has been made in the last decade [28-53] , mainly focusing on the treatment 
of the double real emission^, which resulted in the fully differential calculations of Higgs 
production via gluon fusion [54-57], Drell-Yan [58-62], associated Higgs production with a 
vector boson [63], three jet production from e+e" [64-67] and diphoton production [68]. 

However, further development of methods and new ideas are necessary for efficient 
cancellations of infrared singularities and evaluations of novel two-loop amplitudes in more 
complicated LHC processes. 

With this paper we also take the opportunity to complete the second NNLO applica- 
tion, after the fully differential decay H ^ bb [69], using the method of nonlinear mappings 
to factorize singularities in the double real corrections [70]. The double real contribu- 
tions have often been regarded as the bottleneck of NNLO, and this paper therefore also 
demonstrates the validity of the approach as a method for NNLO corrections in hadronic 
collisions. 

The paper is organized as follows: in section 2 we set up the notation and describe 
the main components of the calculation. In section 3 we provide some detail about the 
treatment of the separation of soft and hard contributions. In section 4 we describe the 
treatment of the double real and the real-virtual components. In section 5 we present 
the way we perform the (non-trivial at NNLO) convolutions for the collinear subtraction 
terms in mass factorization. In section 6 we provide various numerical results both on jet 
rates, pT and rapidity distributions; demonstrate the completely differential nature of our 
calculation and provide typical results for the case in which the Higgs boson decays to two 
photons, including standard experimental cuts on photon momenta and isolation. 

2 Notational setup and conventions 
2.1 Fully differential calculations 

One of the merits of fully differential calculations is the possibility to arrive at theoretical 
predictions for observables in the presence of final state phase-space cuts, like those used in 
experimental analyses, under the precondition that the observable defined is infra-red safe. 
Throughout this article the dependence on such arbitrary phase-space constraints will be 
contained in the jet-function where {p}f denotes the set of final state momenta 

in the lab frame. We will refer to the fully differential cross section as 0"[i7], which we 

variety of methods has been proposed covering the range from fully orthodox to outright heretic. 
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schematically define as 



[J] = Y, dafj{{p}f) 



(2.1) 



/ 



where the sum is over all final states /. 

The usual role of the jet function J is to apply arbitrary final state phase-space 
cuts while ensuring infra-red safety. Here we promote it to a further task, which is to 
keep track of the bin-integrated cross section for any given differential observable with or 
without applying phase-space cuts. This can be achieved simply at the level of Monte Carlo 
integration by passing to J7 not only the set of final state momenta but also the weight of 
the given event. The role of the jet function becomes crucial in all amplitudes that have 
soft and collinear singularities which are regulated by counter terms. In such cases the jet 
function is keeping track of the kinematics of every subtraction term. 

2.2 Hadronic cross section 

We consider the following hadronic process 



where Pi, P2 are the incoming hadrons, H denotes the Higgs boson and X generically 
denotes surplus QCD radiation in the final state. The Higgs boson is assumed to couple 
only to bottom quarks via the SM Yukawa interaction. Assuming the usual factorization, 
the fully exclusive hadronic cross section can be written as 



where the fi{x) denote the bare (unrenormalized) parton distribution functions (PDFs) in 
the 5FS, xi and X2 are the usual Bjorken-x momentum fractions of the partons ii and ^2 

2 

respectively, and r = , where m'jj is the (on-shell) mass of the Higgs-boson and S is the 
square of the total center of mass (CoM) energy of the colliding hadrons. By Ui^i^^HX we 
denote the partonic cross section for the processes 

h{pi) +^2(^2) H{ph) + X{is{ps),u{pi), . . .) , ii,2,3,... e {b,q,g,q,b}. (2.4) 

The PDFs we have inserted in eq.(2.3) are bare and we still have to rewrite them in terms 
of the renormalized PDFs. This step will introduce collinear counter terms that cancel the 
initial state collinear singularities of the partonic cross section, which remain after all real 
and virtual corrections are added together. This cancellation is achieved fully numerically 
in our calculation. We outline the way collinear counter terms can be computed process- 
independently in section 5. 





ap^P2^H+x[J] = / dxidX2 0{xiX2 - T)fi-^{xi)fi2{x2) 



0'iii2^HX [J] 



(2.3) 
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H 



Figure 1. The leading order contribution to bb — > H , see section 2.3. 
2.3 Partonic cross sections 

Expanding the partonic cross section to NNLO in QCD we obtain 



(yij^Hx[J] = Ub 



2 



vr 



(2.5) 



where a sum is impUed over all final state flavors k and / leading to a possible subprocess. 
Here = VbifJ-) and = as{fJ-) are the MS renormalized bottom Yukawa and strong 
couplings, with 5 active flavors. We set the dimensional regularization scale /i to be equal 
to both the renormalization and factorization scales, fipt and fip- Separation of the two 
scales can be easily achieved via the relations given in appendix B. 

The leading order contribution, fig. 1, is denoted by af-^^. At NLO there are two 
separate contributions (see fig. 2): 

• The real {a^j^uk)'- 

Corresponding to the emission of an extra particle k. 

• The virtual [oYj^u): 

Corresponding to the emission and re-absorption of a virtual particle. 

At NNLO there are three separate contributions (see fig. 3): 

• The double real {a^^uj^i): 

Corresponding to the emission of two particles k and /. 

• The real-virtual {crfjY^^j^): 

Corresponding to the emission of one particle k as well as the emission and re- 
absorption of a virtual particle. 

• The double virtual (aYjY^H)'- 

Corresponding to the emission and re- absorption of two virtual particles. 

Real and virtual corrections suffer from infra-red as well as ultra-violet divergences. We 
use conventional dimensional regularization with d = 4 — 2e to regularize such divergences, 
which then appear as poles in e. More specifically ultra-violet divergences present in virtual 
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Figure 2. Some of the diagrams contributing to bb H at NLO. These contributions are denoted 
as real (left) and virtual (right). 



H 




H 



H 



Figure 3. Some of the diagrams contributing to bb ^ H at NNLO. These contributions are 
denoted as double real (left), real- virtual (center) and double virtual (right). 



corrections are absorbed into the renormalized couplings, see e.g. [19, 69] or any textbook 
on QCD. In contrast, infra-red divergences cancel only after summing all real and virtual 
corrections contributing to a given infra-red safe observable. 

Factorized singularities on the unit hypercube may be dealt with using the plus- 
distribution expansion 



-l—ne 



6{x) 



ne ^-^ m\ 

m=0 



(2.6) 



where Vm ix) 



ln"'(x) 



and the plus-distribution is defined through 



V X 



/ dxf{x) 


'aix)' 


-j 


Jo 


X 


+ JO 



(2.7) 



Beyond NLO a factorization of singularities is highly non-trivial. In this work it is achieved 
systematically using the method of nonlinear mappings [70]. 

Care must be taken when dealing with infrared singularities of real emission ampli- 
tudes: the plus-distribution, eq.(2.6), also acts on the jet function, such that cancellations 
happen also at the differential level and are therefore fully local. 



We now give a brief overview of the matrix elements and phase-space measures required 
for the computation of the partonic cross section. Here we assume the amplitudes to be 
color and spin summed. Averaging and phase-space symmetry factors will be explicitly 
factored out. 

i) Purely virtual corrections: 

The purely virtual corrections include the Born, virtual and double virtual contribu- 
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tions and are of the form 

^bi^niJ] = ^J d^2^i\M,-,_j,\^J{p,,p2,PH) , (2.8) 

where S12 = (pi+P2)^ is the partonic CoM energy and N^^i^ = 1/36. The corresponding 
phase-space volume element is trivially given by d^2^i = 27r(5(si2— m|^), constraining 
si2 = m'jj = p'jj. Regarding the computation of amplitudes we refer the reader to [69] 
where the full virtual matrix-elements can be found. Explicit expressions for these 
contributions are given in section 4.4. 

ii) Single real emissions: 

The single real emissions include real and real- virtual corrections and are of the form 

N- ■ f 

ai^i2-,Hi:i[J] = TT^ / d^2^2\Mi^i^^Hi3\'^JiPl,P2,P3,PH)- (2.9) 
^Si2 J 

Further details are given in section 4.1. 

iii) Double real emissions: 
These are of the form 

ai^i^^HiauiJ] = j d^2^3\Mi^i2^Hi3ii\'^JiPi,P2,P3,P'i,PH)- (2.10) 

Further details are given in section 4.3. 

3 Separation of soft and hard 

Since in the soft limit of a 2 —t- 1 process the produced particle is at rest in the partonic 
center of mass frame, there is no difference between the soft piece of a fully differential 
partonic cross section and that of a fully inclusive partonic cross section. It is therefore 
very convenient to isolate the soft contribution (crs) to the partonic cross section (a) from 
the hard one {an), i-e. 

a = as + crH. (3.1) 

This allows for a fully analytic treatment of as, while an must, as far as external kinematics 
are concerned, be treated numerically. Let us introduce the variable 

2 

z = —^. 3.2 

Sl2 

Then the soft limit of all real emission amplitudes corresponds to z — )• 1, which identifies 
the production threshold. Given that infrared singularities are of logarithmic nature, the 
divergence at z = 1 can be exposed as follows 

a{z)[J] = 6{l - z)av{e)[J]U=, + J] '^i^^ , (3.3) 

where ay denotes the purely virtual correction, while a^^\z,e) denotes real corrections 
collectively (at NNLO this includes both real- virtual as well as double real corrections). 
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Separation into soft and hard parts can now be achieved by adding and subtracting 
the soft hmit from the second term in the above, yielding 



such that (Th is integrable in the range z G [t, !]• Of course this decomposition of the 
partonic cross section into its soft and hard components is not unique: one could use 
any other subtraction term with the correct limit, thereby including, for example, the 
luminosity function. Our choice, however, has the nice property that the soft part as can 
be expanded purely in terms of 6- and plus-distributions via eq.(2.6). 



TS[J] = Co5(l - Z) J\,=l + V CnVnil " z) J\ 



n=0 



2 = 1 



~ ~(n) 

Thereby all threshold divergences between ay and cr]^ are canceled analytically, leaving 
only a finite threshold contribution. Furthermore this framework provides a natural way 
to incorporate threshold re-summation in fully differential calculations. 

4 Details on the calculation 
4.1 The single real 

The single real partonic cross section may be expressed as 

[^] = y" (v) £ / d^^\M!!\'J{Pi,P2,P3,PH) . (4.1) 

We define d^2 to be the conventional phase-space volume d^2 up to some renormalisation 
constants. Here we have to consider 6 separate channels: 

'r^lMi^nnl' if(i,i)G{(6,6),(6,6)}; 



bb^gW 

H2-ly,., \Kg^bH? if ihJ) e {(^9), (6,5)} 

tR 



Kb-.bHV if(i,j)G{(5,&),(5,6)}. 



, 2-{2-2e)-3-8 \^'^^gb-^bH 

The corresponding amplitudes may all be found in [69] . A convenient phase space parametriza- 
tion is given by 
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where A G [0, 1], with the Lorentz invariants taking the simple form 



si3 = (Pi - Psf = -si2(l - z)X , 

S23 = iP2 - Psf = -si2(l - ^)(1 - A) . (4.3) 

Note that the singularities of S13 and S23 are factorized in A, (1 — A) and (1 — z) which 
allows for a simple subtraction of the poles using eq.(2.6). This also allows us to identify 

""hh Vf\ = (1 _ ^y+2e ■ (4.4) 
The calculation of the hard part then trivially follows from eq.(3.4). 
4.2 The real-virtual 

The real-virtual partonic cross section may be expressed as 

'^^'^ [^] = (2^^)'(v)'£ / d^22Re{M,fM,f}jip,,p2,ps,PH), 

where we have taken the liberty to define d^2 to equal eq.(4.2) up to some renormalization 
constants. Then 



iV,j2Re {m^^M^*} 



2.(2-L).3.8 if(^,j)e{(6,<7),(6,5)}; 

k 2.(2-2e).3.8 )1 if(^,j)e{(5,&),(9,W- 



The real-virtual amplitude can be obtained from the corresponding one from the decay 
process H ^ bb published in [69] by crossing particles to the initial state. The box integrals 
we encounter in this amplitude are entirely expressible in terms of Gauss' hypergeometric 
function 2-^1(1, — 1 — z) where z can be in any of the three sets Sjim, Sinv and Sni '■ 

q - / ~^13 -S23 \ _ / -S12 -S12 

•J fine — ) ) ( } Jinv — ) i 

[ Si2 Si2 J [ Si3 S23 

S 1 = 1-^ -si2mjj -sumjj -S2?.m]j \ ^^^^ 

I S23 ' si3 ' S13S23 ' S23S12 ' S12S13 J ' 

When attempting a direct subtraction of the singularities created by the real emission, 
the points of subtraction overlap with singular points of the hypergeometric functions in 
the box integrals. It was found in [69] that one can apply transformations on the argument 
of the functions to circumvent this difficulty. Since here we are no longer in the euclidean 
regime of this amplitude, the required transformations are different than in [69]. Analyzing 
integral representations, we find that we have to apply the following identities: 

• If z G Sfine the soft-collinear limits are well defined. 
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li z G Sni we apply 

2^1(0, b; c; z) !-)• (1 — z)~^2Fi [ c — a,b;c; 



z - 1 



• If z G Sinv we employ the argument inversion, 

T{b-a)T{c)2Fi{a,a-c + l;a-b+l;l) 

^^^("''''^'^^^ r(5)r(c-a)(-.)" 

r(a-b)r(c)2Fi(fe,b-c + l;fe-a + l;^) 

r (a) r (c - 6) (-z)^' ■ ^ ■ ^ 

After these transformations are applied, the singularities corresponding to the real emission 
are factorized in A, (1 — A) and {1 — z). The soft singularity structure of the real-virtual 
may then be extracted as 

4 ~{m)RV r 

In the soft limit only the m = 2, 4 coefficients survive and the integration over A can be 
done analytically. The explicit expressions for the soft limit can be found in appendix A. 

The computation of the hard part then follows from eq.(3.4). While the structure 
is more complicated than in the case of the single real, a direct subtraction via eq.(2.6) 
can still be achieved in a straightforward manner. In order to obtain the final Laurent 
expansion in e we employ the library HypExp [71] to expand the hypergeometric functions 
in terms of polylogarithms. 



4.3 The double real 

The double real partonic cross section can be written as 

where d^s is equal to the conventional three-particle phase space element up to renor- 
malization constants. Using the discrete symmetries of the squared amplitudes we are able 
to considerably reduce the number of independent channels, which one has to implement 
separately. These symmetries are due to the charge invariance of all the bb ^ H double 
real amplitudes (exchanging g -f-)- g or 6 o 6 leaves the amplitudes invariant). This leaves 
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us with the following list of channels 



RR\2 



til -yi [ti l-^fif,— s-oof/l + 










if 




e{{b,b),{b, b)}; 




if 


ihj) 


e {{q,q), {q,q)} ; 


1 ]\dr,h ^nh U 1 ^ 

ozoz l-^^-^gO — rqOJri \ 


if 


i. 


e Uci,b), (Q,b), fa, 5), fo,6)l 


2^\Mbq-^bqH\ 


if 




G{(6,(7),(6,g),(&,g),(6,g)} 


2-(2-2e)-3-8 


if 


ihj) 


e{(6,5),(^,<7)}; 


1 \nj |2 


if 


(hj) 


(^{{9,b),{g,b)}; 


1 Uf - |2 


if 


(iJ) 


= {9,9) ; 


27?37 2! ^''^'-s^b^'-fi' 


if 




G {(6,6), (6,6)}. 



(4.8) 

By crossing partons from the initial to the final state, we can obtain all of the above from 
the three amplitudes 



IM, 



IM, 



bbbbH^O 



P and IM 



bbqqH^O I 



published in [69]. In order to deal with the intricate singularities, their factorization and 
subtraction, we refer the reader to the methods developed in [70], which we have imple- 
mented faithfully. As in the single real emissions, the double soft singularity occurs at the 
threshold. Its structure may be identified as 



RiJr^l "bb L^J 



(l-z) 



l+4e 



(4.9) 



4.4 The soft 

Let us expand as in the strong coupling 



as[J] =)3-{6{l-z) + ^As,NLO + (- 



A 



S,NNLO 



+ 0{al)] J\ 



z=l 



where 



The NLO correction A 



(4.10) 
(4.11) 



S,NLO 



may be expressed as 



As,NLO 



8 16 

+ -Vo{l-z)lH + -Viil-z) + Oie) 



5(1- z) 



while the NNLO correction As^nnlo can be expanded as follows 



(4.12) 



(4.13) 
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with the only non-vanishing contributions [72] 

a(0) { { 10^ 2 \ 64^,2 / 17 8^ 250 \, 

^s'nnlo =[\jj~ ^"^2 + 3C3 \nf- —C2IH + ( ~y + 3^2 + ^Cs 1 Ih 



211 58^ 26^ 17 ^V,, 

H ^ C2 Cs C4 M (1 - ^) 

18 9 ^ 3 ^ 6 ^ V ^ 

/ /56 20, 2, 2 8^ \ ,2 /70 164^ \ , 

212 20^ 638^ \^ ,^ 
"^ + y^' + "9"^V °^ 

+ ((^^.-f)n, + f + lf/.^-4..-72C.)p.(l-.) 

+ {-A + ^^Uf + ^/h) V, (1 - z) + (1 - ^) , (4.14) 



(-1) / / 1 2 ^ \ 64 _ 43 23 . 125 



^In'nLO = ^ (^^ + 9C2J -/ + ^C2lH + ^ - y C2 - ^CsJ <^ (1 - Z) 

/lO 16, 35 82 \ _ 
I + "^(ff) Pi (1 - .) - ^1)2 (1 - .) (4.15) 



and 

^^NNLO = (t " l^^f ~ f ^0 5 (1 - ^) + (9 - ^n;) Po (1 - z) 

+ ^Pi(l-z). (4.16) 
Here n/ is the number of hght flavors, Cn are the usual Riemann zeta values and 

Ih = log (^) . (4.17) 

The explicit soft limits of the real-virtual and double real pieces that are included in its' 
can be found separately, and with their full color-dependence, in appendix A. 

5 Collinear factorization 

Parton distribution functions are renormalized to absorb initial state collinear singularities 
via 

~fi{z,^l) = [Vi,{^x)®fj){z), (5.1) 

where /i is the factorization scale and fi are the bare parton densities. In the following 
discussion summation over indices will always be assumed unless explicitly stated. We will 
also need the convolution integral, which is defined as 

U®9){z)= I dxdyf{x)g{y)5{z-xy). (5.2) 
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The kernel Fj,- is defined in the MS scheme by 



r,,(z) = 6.,5{l - ^) + ( f ) Tf^{z) + )' rg)(^) + 0{al) , (5.3) 



where the coefficients of the expansion in the strong couphng involve the Altarelli-Parisi 
splitting functions P^^. Specifically, 

rg)(z) = -^, (5.4) 

rlf (^) = -[^~^2[{P^® Pk,) (-) + PoPU^)] I , (5.5) 

with /3o = X ~ define the inverse of the kernel Fjj as 

A,,(.) = ^AS;)(.)(^)%0(a3), (5.6) 

n=0 

such that it satisfies the condition (Fjfe (g) A^j) (z) = 5ij5{l — z). Solving for the coefficients 
yields 

(5.7) 
(5.8) 

r«^r«)(.) 

= ^ + 2^ [(^^i ^ P'ki) (-) - P^PfM ■ (5.9) 
The strong coupling expansion of the bare PDFs then reads 

M-) = E/f^(-)(v)"+^K)' (5-10) 

n=0 

with 

/f^ = a1;) /, . (5.11) 

In evaluating the collinear counter terms we encounter convolutions of the type (/ (g) A) (x) , 
where the function / is regular and A(a:) can in general be written as 

A(x) = a5{l -x) + Y, bnDn{x) + C{x). (5.12) 

n 

Expressing the convolution as a single integral we obtain 

(/ ® A)(x) = £ ^ / (^) |«'5(i - y) + E ^"^"(y) + ^(y) I • (5-13) 



aS?(.) 


= 5ij5{l - z) , 




Ag^(.) 


= = 




Ag^(.) 


= -rg'(.) + 
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Care has to be taken with convolutions over T>n- Since the integration does not start at 
zero, a boundary term must be included 



log(l - 



f{x) + dy log(l - y) 



m 



1 

y 



/(f)-/(x) 



(5.14) 



m + 1 



X 



Because of the downward sloping shape of all parton distribution functions, a quadratic 
remapping of the integration variable y was found to optimize the convergence behavior, 
i.e. we parametrized the integral like 



with z uniformly distributed between and 1. 

In our code, this integration is carried out numerically. The integration is one- 
dimensional, which makes a simple deterministic trapezium integration with about 50.000 
points the simplest option. The result of the integration is accurate to at least 5 digits, 
which is usually below the precision of the Monte Carlo integration. The precision of the 
integration can be arbitrarily increased by increasing the number of points used. For ev- 
ery bare PDF used, we construct a one-dimensional grid in the Bjorken-x variable and 
interpolate from it during runtime. 

An alternative to constructing a grid is to perform the integration numerically along 
with the phase space ones, thereby increasing the dimensionality of the Monte Carlo inte- 
gration by one (or by two in the case of double NLO kernels convoluted with the Born). 
We have implemented this as well and found that it yields the same results as the grid 
approach. 

This procedure allows us to expand the (singular) bare PDFs via eq.(5.10) order by 
order in the dimensional regulator e and substitute them directly in eq.(2.3). The singular- 
ities in the resulting convolutions, appearing as poles in the e-expansion, cancel the initial 
state collinear singularities of the partonic cross section. This cancellation is achieved nu- 
merically in our calculation and can be observed bin by bin in e.g. the rapidity distribution 
of the Higgs boson. One can achieve this cancellation in each initial state channel sepa- 
rately, at the cost of separating the convolution integrals depending on the initial state 
parton in the convolution, i.e. by not performing the implicit j-summation in eq.(5.11). 

It is worth pointing out that the procedure described here is entirely generic, i.e. it 
provides the collinear counter terms for any NNLO process numerically. Moreover, we 
thereby circumvent the usual insertion of eq.(5.1) in the equivalent of eq.(2.3) for renor- 
malized quantities and the resulting cumbersome and process specific analytic treatment 
of the convolutions. 

6 Numerical results 

We have performed a number of tests to ensure that our results are consistent with each 
other and with results available in the literature: 
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• We have implemented the entire calculation in two different computer codes, one in 
Fortran and one in C++, and all results agree within their respective Monte Carlo 
errors, both inclusively and differentially. 

• The coefficients of all poles in the e-expansion of all cross sections cancel both inclu- 
sively and differentially for the entire process and also for all individual initial state 
channels. 

• The inclusive cross section agrees with the one available in [19] and from ihixs [73] 
and so does the inclusive cross section per initial state channel. This is the first 
independent check of the inclusive cross section published in [19] and adopted in [73]. 

• The soft limit of both real- virtual and double real contributions were computed both 
numerically (as a limiting case of the generic matrix elements) and analytically. More- 
over the integrated double real contributions were found to agree with an analytic 
computation provided by [72] . 

• The subtraction process for every double real integral was implemented in two dif- 
ferent ways and were found in complete agreement. 




Figure 4. The Higgs rapidity distribution for mn = 125 GeV at the 8 TeV LHC. The bands 
describe the uncertainty due to factorization scale 

We present results for the LHC with a center of mass energy of 8 TeV. We fix the 
mass of the Higgs boson at 125 GeV. We have used the MSTW2008 (68%CL) PDFs for aU 
results presented here. The value of Og at mz that we use is the best-fit value of the PDF 
set at the corresponding order. We use = mn as the central renormalization scale. The 
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Figure 5. The Higgs transverse moraentum distribution for tjih — 125 GeV at the 8 TeV LHC. 

value of Us used is run from mz to through NNLO in QCD. The mass of the bottom 
quarks is set to zero in ah matrix elements, consistently with the 5FS choice. The bottom 
Yukawa coupling, however, depends on the mass of the bottom. The Yukawa coupling at 

is obtained from the Yukawa coupling at /U* = 10 GeV, using mh{n*) = 3.63 GeV. 

We do not vary fiR in what follows, since the scale dependence of the total cross 
section has been found to be very mild. We have also checked that the //^-dependence of 
differential distributions is very small. 

Previous studies have shown that the inclusive cross section is very sensitive to the 
choice of factorization scale. Arguments related to the validity of the 5FS approximation 
with respect to the collinearity of final state 6-quarks, as well as to the matching to the 4FS 
calculation or to the need for a smoother perturbative expansion, point to factorization 
scales that are much lower than the Higgs boson mass. We adopt the choice fip = as 
a central scale and vary it in the range [^^, to estimate the related uncertainty. 

All Monte Carlo integrations was performed with the Cuba [74] implementation of the 
Vegas algorithm. 

The rapidity distribution of the Higgs boson is shown at fig. 4. As expected, the 
perturbative expansion is converging smoothly for this choice of central fip and the NNLO 
uncertainty band is entirely engulfed by the NLO one. 

The transverse momentum distribution for the Higgs boson is shown in fig. 5. This 
observable starts at NLO in QCD in the 5FS, and the fixed order prediction fails, as usual, 
to describe the very low pT spectrum due to the related large logarithms. At the large px 
range we see that the NNLO calculation leads to a harder spectrum than the NLO one 
and the NLO scale uncertainty fails to capture this feature. This implies that great care 
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should be taken when relying on NLO predictions for observables that are highly exclusive 
in the transverse momentum of the Higgs boson. 




Figure 6. Differential distribution in rapidity and transverse momentum of the Higgs boson for 
tuh = 125 GeV at the 8 TeV LHC. 

The differential distribution in both the rapidity and the pT of the Higgs is shown in 
fig. 6, both in a three-dimensional lego plot and in a density plot. We see that the bulk of 
the events are produced centrally (with \y\ < 2.5) and at relatively low px { 35 — 50GeV). 




Figure 7. The cumulative distribution of the Higgs pr for iriH — 125 GeV at the 8 TeV LHC. 

In fig. 7 we show the cumulative distribution of the Higgs transverse momentum. This 
observable is equivalent to the cross section in the presence of a jet veto at NLO, but only 
related indirectly at NNLO. In fig. 8 we present the cross section in the presence of a jet 
veto. We see again that the perturbative description for high pT cut-offs is satisfactory 
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Figure 8. The cross section in the presence of a jet veto (the 0-jet rate) for mn = 125 GeV at the 
8 TeV LHC. 



(despite the discrepancy in high px between NLO and NNLO, which is, in absolute terms, 
unimportant), while for cut-offs lower than 20 GeV the NLO description does not coincide 
with the NNLO one. The vanishing of the uncertainty around 15 GeV (which in the case 
of the jet veto is taking place at a slightly lower p^-veto value) is a feature reminiscent of 
a similar situation in Higgs production via gluon fusion [75]. The fixed order prediction in 
this region is very stable under varying the factorization scale, and any residual uncertainty 
in quantities like the acceptance in the presence of a veto is driven by the uncertainty in the 
total cross section. Various approaches to assign a larger uncertainty to similar observables 
involving re-summation exist, see for example [76]. 

An important observable in bb ^ H is the cross section for zero, one and two jets. We 
use the anti-Zc-r algorithm [77] for jet clustering^ with a cone in the y — (j) plane of radius 
R = 0.4. We show in fig. 10 the jet rates as a function of the jet used to define 

them. Here we do not distinguish between 6-jets and light jets. We find the jet rates for 
p-max _ 20GeV to be in agreement with those published in [24]. 

A wealth of information can be derived from examining the contribution of the differ- 
ent initial state channels to differential distributions. The six initial state channels that 
contribute to our NNLO calculation have singularities in various collinear regions that are 
canceled against the collinear counter terms from mass factorization. In order to make the 
cross section per channel finite one has to use collinear counter terms that include rj™^ 
kernels involving only the initial state partons of the channel considered. Since we calculate 

^At this order in perturbation theory, the anti-fcr, the kr and the Cambridge-Aachen algorithms are 
completely equivalent. 
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Figure 9. The average pt (left) and rapidity (right) of the Higgs as a function of fip/mn for 
niH = 125 GeV at the 8 TeV LHC. 




Figure 10. The 0-, 1- and 2-jet rate as a function of the pt used in the jet definition for niH — 125 
GeV at the 8 TeV LHC. 

the collinear counter terms numerically this modification was relatively easy to achieve. 

Initial state channel contributions to differential distributions have a strong dependence 
on the factorization scale, as do initial state channel contributions to the inclusive cross 
section. In fig. 11 we see the contributions to the Higgs boson px distribution from each 
channel, for various factorization scales ranging from m/^/16 to 2m//. 

Within the 5FS, the factorization scale regularizes the collinear singularities which in 
the 4FS are regularized by the bottom mass. At NNLO, three initial state channels, bb, bg 
and gg share common collinear configurations whose leading logarithms cancel each other 
in different bins of the Higgs pT distribution. In the zero px bin, in particular, squared log- 
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Figure 11. The distribution of the Higgs pT per initial state channel for mn = 125 GeV at the 



TeV LHC, with fip^"^. 



rriH niH mjj 



arithms from the double collinear limit of the gg channel cancel against the single collinear 
limit of the bg channel and the born contribution of the bb channel. Moreover at NNLO one 
also sees sub-leading (single) logarithms canceling each other between the single collinear 
configurations of the gg channel and the regular contributions to the bg channel, a can- 
cellation that appears in non-zero pT bins as well. The magnitude of those logarithmic 
cancellations is regulated by the value of the factorization scale. The factorization scale 
dependence is an artifact of the truncation of the perturbative series, so one would naively 
choose the scale in a way that minimizes the cross-channel logarithmic cancellations. How- 
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Figure 12. The distribution of the Higgs absolute rapidity, \y\ per initial state channel for ttih = 
125 GeV at the 8 TeV LHC, with fip = '^,mH,2mH. 

ever, choosing the scale too small reduces the regime where the logarithms are re-summed 
in the PDFs, destabilizing the perturbative expansion. Ideally one should choose the scale 
in the region where the collinear approximation implicit in the 5FS is still reasonable, which 
is at m._f//4 or lower. Corroborative evidence for such a choice comes from the behavior 
of the average transverse momentum and of the average rapidity of the Higgs boson as a 
function of the factorization scale choice, shown in fig. 9. 

These features are also seen in the rapidity distribution of the Higgs boson per initial 
state channel, shown in fig. 12 for various values of ^p- There it is clearly seen that a scale 
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Figure 13. The average px of the two photons and the Y* distribution for bb—>-H + X—>-jj + X 
for niH = 125 GeV at the 8 TeV LHC, in the presence of cuts described in the text. 

like fiF = f^n/^ eliminates the cross-channel cancelations but a lower scale iip = 
leads to a reduced, 6g-dominated prediction. 

We turn now to more exclusive observables. In large tan /3 models where the Higgs 
boson production gets significant contribution from the bottom quark annihilation process, 
one would like to examine differential distributions involving decay products of the Higgs 
boson, with cuts necessary in the experimental analyses. We focus here, for demonstration 
purposes, on the case where the Higgs boson decays to two photons. In such an analysis 
the minimal cuts used by CMS and ATLAS include: 

• A cut on the pT of the leading photon: px-^i > 40GeV. 

• A cut on the pT of the trailing photon: pT-2 > 25GeV. 

• A cut on the rapidity of both photons: |yi,2| < 2.4. 

• An isolation cut on photons: no jet is allowed in a cone of radius 0.4 around any of 
the two photons if it is px > 15GeV. 

We treat the Higgs boson in the zero width approximation in this article. We defer a more 
realistic treatment of the Higgs propagator to future work. 

Within this setup we show in fig. 13 the distribution of the average transverse momen- 
tum of the two photons and the distribution of the absolute of the difference in pseudo- 
rapidity between the two photons, Y* = ^\yi — 2/2 1- 

7 Conclusions 

We have presented the fully differential NNLO calculation of bb ^ H, a process of prime 
phenomenological importance for the LHC in all models with enhanced bottom Yukawa 
couplings. This is the first independent cross-check of the inclusive NNLO calculation 
performed in [19]. We have presented a variety of differential distributions for Higgs pro- 
duction that can only be obtained with a fully differential calculation and are useful for 
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assessing the quality of the perturbative expansion and the level under which several fea- 
tures are under control at a fully differential level. We have also presented predictions for 
fully exclusive observables for the 66 — >• — t- 77 process in the presence of tight cuts on the 
final state photons including isolation cuts, demonstrating that our calculation can fully 
simulate any experimental setup at the partonic level. 

This is the second application of our approach to treat real emission singular amplitudes 
at NNLO [70] . It is the first application for the more complicated case of a hadron collider 
process. We find the approach particularly beneficial, both in terms of automatization and 
in terms of performance of the resulting numerical code. We find that the improvement in 
performance compared to the sector decomposition approach is significant. We intend to 
release the computer code in the near future and we defer for then any detailed comments 
on performance issues. 

A study of significantly wider scope, including the production via gluon fusion in 
models with enhanced bottom Yukawa couplings, as well as the decay of Higgs to bottom 
quarks or tau leptons would vastly benefit the experimental searches. We defer such a 
study for a future publication. 
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A Threshold behaviour 
A.l Double real 

The z dependence of the double real soft contribution factorizes completely 

af« = 4CF^e-2'«^ (1 - Agi, 
where lu is defined in eq.(4.17), B in eq.(4.11) and 

/ 1 _2 5 _i 7 7 ^ / 41 35^ 31 \ \ 
+ +72^ + 54 - 24^^+Vl^-72^^-36^V7"' 

/3 „3 11 _2 / 67 11 \ _i 101 77. 67^ 
"U' ^48' +Vn4-T^V' +l08"48^^"¥^^ 
/607 469^ 341 19^ \ \^ 
+ (324 - 144^^-T^^^-32^VV'^" + ''^^^- 
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A. 2 Real- virtual 

We decompose the real-virtual soft contribution as 



where only A^y and are non vanishing and given by 
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B Scale separation 

The rcnormalization and factorization scales, /i/j and can be conveniently separated 
by first setting ^ = and then applying the following relations 
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